function [ T, area, Ix, Iy, xb, yb ] = calcMomentOfAreaTensor(contour)
% function [ T, area, Ix, Iy, xb, yb, C, c_H, c_var, e_var, r0 ] = calcMomentOfAreaTensor(contour)
N = size(contour,1);
%% area, center-of-mass
accArea = 0;
accIx = 0;
accIy = 0;

for k = 0:N-1
    l = mod(k+1,N);
    xi = contour(k+1,1); xj = contour(l+1,1);
    yi = contour(k+1,2); yj = contour(l+1,2);
    
    accArea = accArea + xi*yj - xj*yi;
    accIx = accIx + (xi*yj - xj*yi)*(xi + xj);
    accIy = accIy + (xi*yj - xj*yi)*(yi + yj);
end
area = 1/2 * accArea;
Ix = 1/6 * accIx;
Iy = 1/6 * accIy;

xb = Ix/area;
yb = Iy/area;
%% moment-of-area tensor
accT11 = 0;
accT22 = 0;
accT12 = 0;
for k = 0:N-1
    l = mod(k+1,N);
    xi = contour(k+1,1)-xb; xj = contour(l+1,1)-xb;
    yi = contour(k+1,2)-yb; yj = contour(l+1,2)-yb;
    
    accT11 = accT11 + (xi*yj - xj*yi)*(yi^2 + yi*yj + yj^2);
    accT22 = accT22 + (xi*yj - xj*yi)*(xi^2 + xi*xj + xj^2);
    accT12 = accT12 + (xi*yj - xj*yi)*(2*xi*yi + xi*yj + xj*yi + 2*xj*yj);
end
T11 = 1/12 * accT11;
T22 = 1/12 * accT22;
T12 = 1/24 * accT12;
T = [T11 T12; T12 T22];
% %% 2. order moments
% x = contour(:,1)-xb;
% y = contour(:,2)-yb;
% 
% c_xx = 1/N*sum(x.^2);
% c_xy = 1/N*sum(x.*y);
% c_yy = 1/N*sum(y.^2);
% C = [c_xx c_xy; c_xy c_yy];
% 
% % c_rpax = (c_yy + c_xx - sqrt( (c_yy + c_xx).^2 - 4*(c_xx.*c_yy - c_xy.^2) )) ./ (c_yy + c_xx + sqrt( (c_yy + c_xx).^2 - 4*(c_xx.*c_yy - c_xy.^2) ));
% %% modified haralick's circularity
% r = sqrt(x.^2 + y.^2);
% c_H = mean(r) / std(r);
% %% circular variance
% c_var = 1/(mean(r)^2) * var(r,1);
% %% elliptic variance
% C_inv = inv(C);
% r_ell = sqrt( C_inv(1,1)*x.^2 + (C_inv(1,2)+C_inv(2,1))*x.*y + C_inv(2,2)*y.^2 );
% e_var = 1/mean(r_ell) * var(r_ell,1);
% %% front radius
% contour = sortrows(contour);
% num = 6;
% [~,~,r0] = circfit(contour(1:num,1),contour(1:num,2));
% % contourCircle = sortrows(contour(1:num,:),2);
% % m1 = -(contourCircle(2,1)-contourCircle(1,1))/(contourCircle(2,2)-contourCircle(1,2));
% % n1 = (contourCircle(1,2)+contourCircle(2,2))/2 - m1*(contourCircle(1,1)+contourCircle(2,1))/2;
% % m2 = -(contourCircle(4,1)-contourCircle(3,1))/(contourCircle(4,2)-contourCircle(3,2));
% % n2 = (contourCircle(3,2)+contourCircle(4,2))/2 - m2*(contourCircle(3,1)+contourCircle(4,1))/2;
% % xM = (n2-n1)/(m1-m2);
% % yM = m1*xM + n1;
% % x = contourCircle(:,1) - xM;
% % y = contourCircle(:,2) - yM;
% % r = sqrt(x.^2 + y.^2);
% % r0 = mean(r);
% % opts = optimset('MaxFunEvals',50000, 'MaxIter',10000, 'TolFun',1e-20, 'TolX',1e-5, 'Display','off');
% % para0 = [5,5,10];
% % xData = contour(1:num,1); xDataShift = min(xData);
% % yData = contour(1:num,2); yDataShift = min(yData);
% % dataOLS = @(para) sum((sqrt((xData-xDataShift-para(1)).^2 + (yData-yDataShift-para(2)).^2) - para(3)).^2); % Ordinary Least Squares cost function
% % dataParaOpt = fminsearch(dataOLS, para0, opts);
% % % x0 = dataParaOpt(1)+xDataShift;
% % % y0 = dataParaOpt(2)+yDataShift;
% % r0 = dataParaOpt(3);
% % % phiFit = linspace(0,2*pi);
% % % xFit = r0*cos(phiFit) + x0;
% % % yFit = r0*sin(phiFit) + y0;
end